/**********************************************************************
*math.h - definitions and declarations for my math library
*
*Author : Zhuowen Tu
*
*Purpose:
*       This file contains constant definitions and external subroutine
*       declarations for the math subroutine library.
*
*       [Public]
*
***********************************************************************/
#ifndef _MC_PUBLIC_H
#define _MC_PUBLIC_H

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <io.h>
#include <iostream>


#include "mcPoints.h"
#include "mcPoints.h"
#include "mcPoints.h"

#define EBSHLON			1e-10								//definition of minim value
#define	INFINITE_D		1e+99								//definition of maxim value
#define	INCH_TO_METER	0.3048

#ifndef PI
#define	PI				3.1415926535897931
#endif

#define INTENSITY_DELTA		6.0
#define INTENSITY_LOW		35.0
#define INTENSITY_HIGH		255.0

#define DIS_SMALL		0.05

typedef class DpMcPoint3D Coord3D;
typedef unsigned char uchar;
typedef McPoint2D*	PMcPoint2D;

inline double round(const double a)
{
	double b,c;

	b = floor(a);
	c = fmod(a,1.0);
	if (c>=0.5)
		b = b+1;
	return b;
}

//get the not less integer of a floating value
inline int	NotLessInt(double Value)
{
	int Rtn;

	Rtn = (int)Value;
	if (Rtn < Value)
		Rtn++;
	return Rtn;
}

//crossing situation of two lines
#define			CROSSING		1
#define			DONT_INTERSERT	0
#define			PARALLE         -1
#define			ONE_POINT		-2

//calculation of distance between two points
double	Distance(double x1, double y1, double x2, double y2);
double	Distance2(double x1, double y1, double x2, double y2);
inline double	Distance(const McPoint2D xy1, const McPoint2D xy2) {return Distance(xy1.X(), xy1.Y(), xy2.X(), xy2.Y());};
inline double	Distance(McImagePoint xy1, McImagePoint xy2) {return Distance(xy1.X(), xy1.Y(), xy2.X(), xy2.Y());};
double  Distance(double x1, double y1, double z1, double x2, double y2, double z2);
inline double   Distance(McPoint3D xyz1, McPoint3D xyz2) {return Distance(xyz1.x, xyz1.y, xyz1.z, xyz2.x, xyz2.y, xyz2.z);};

//calculation of two lines
int		LineIntersection(McPoint2D *ppt2dCrossing, McPoint2D xy11, McPoint2D xy12,
						 McPoint2D xy21, McPoint2D xy22);

//calculation of perpendicular distance between a point and a line
double  PerDistance(const McPoint2D &pt2d, const McPoint2D &ptLnBg, const McPoint2D &ptLnEd);
double  PerDistance2(const McPoint2D &pt2d, const McPoint2D &ptLnBg, const McPoint2D &ptLnEd);
bool	PointOnLeftSide(const McPoint2D &pt2d, const McPoint2D &ptLnBg, const McPoint2D &ptLnEd);
bool	PointInLine(const McPoint2D &pt, const McPoint2D &pt11, const McPoint2D &pt12);

//calculation of the angle between x axis and line
double	xysa(double x1, double y1, double x2, double y2);

//calculation of the angle between inputed two lines
double	AngleOfTwoLines(McPoint2D xy11, McPoint2D xy12, McPoint2D xy21, McPoint2D xy22);
//calculation of area
double	area(double *x, double *y, int cnt);

//cut blanks in string
void	CutBlank(char * lpszDest, char * lpszSour);
void	CutEndBlank(char * lpszDest, char * lpszSour);
void	CutHeadBlank(char * lpszDest, char * lpszSour);
void	GetFirstWord(char * lpszDest, char * lpszSour, char cSeprator, bool bChangeSour=false);
void	GetLastWord(char * lpszDest, char * lpszSour, char cSeprator, bool bChangeSour=false);

//point stays on input line or not
bool	PointOnLine(McPoint2D point, McPoint2D ptVertex1, McPoint2D ptVertex2);

//point lies inside or ouside the given polygon
#define OUTSIDE			0    //outside
#define INSIDE_CW		1    //inside of a clockwise polygon
#define INSIDE_CCW		2	 //inside of a counterclockwise polygon
int		VeriPoint(McPoint2D *pVertices, int iVrtcCt, double x, double y);
bool	IsLineCrossPlgn(const McPoint2D ptFirst, const McPoint2D ptEnd,
						McPoint2D* pPts, const int iPtCnt);

inline	double	InchToMeter(double a){return a*12*0.0254;};

void	DoubleToString(char *buf, double value);
double	ComputeNormalCurvature(const McPoint2D& pt1, const McPoint2D& pt, const McPoint2D& pt2, double* pcurvature=NULL);
void	PolarCoordinatesLines(double &rou, double &theta, const McPoint2D& pt1, const McPoint2D& pt2);


#define min(a, b)  (((a) < (b)) ? (a) : (b))
#define max(a, b)  (((a) > (b)) ? (a) : (b))

template <class Tp>
int Sign(Tp s)
{
	if (s>0) return 1;
	if (s<0) return -1;
	return 0;
}

// arctg function
inline double Atan2(const double y, const double x)
{
	double theta=atan2(y,x);
	if (theta<0)
		theta += 2*PI;
	return theta;
}

// perpendicular angle
// x 0~2*PI
inline double PerPenAngle(const double x)
{
	double theta = PI/2+x;
	if (theta>=2*PI)
		theta -= 2*PI;
	return theta;
}

double AngleOPI(const double alfa);
double AngleO2PI(const double alfa);


template <class Tp>
Tp ABS(const Tp x)
{
	if (x>0) return x;
	else	 return -x;
}


inline double factorial(const int i)
{
	int k;
	double d=1.0;

	if (i>0)
	{
		for (k=1; k<=i; k++)
			d = d*((double)i);
	}
	return d;
}

inline double choose_i_k(const int i, const int k)
{
	double a1,a2,a3;

	if (k==1 || k==i-1)
		return (double)i;
	else if (k==2)
	{
		return ((double)i)*(i-1)/2.0;
	}

	a1 = factorial(i);
	a2 = factorial(k);
	a3 = factorial(i-k);

	return a1/a2/a3;
}

inline int IntClose(double x)
{
	double y1 = (double)ceil(x);
	double y2 = (double)floor(x);
	if (fabs(y1-x)<fabs(y2-x))
		return (int)y1;
	else
		return (int)y2;
}

inline double IntToDouble(int i)
{
	double y1=floor((double)i);

	if (IntClose(y1) == i)
		return y1;
	else
	{
		y1 = y1+1.0;
		if (IntClose(y1) == i)
			return y1;
		else
		{
			y1 = y1-2.0;
			if (IntClose(y1) == i)
				return y1;
		}
	}
	return y1;
}


inline double ProbFromEnergy(const double energy)
{	if (energy>700)	return 0.0;
	else			return exp(-energy);
};


template<class Tp>
void SwapInt(Tp &x1, Tp &x2)
{
	Tp temp;
	temp = x1;
	x1 = x2;
	x2 = temp;
}

inline double DistanceAngles(const double angle1, const double angle2)
{
	return min(AngleO2PI(angle1-angle2), AngleO2PI(angle2-angle1));
}

inline double AverageAngles(const double angle1, const double angle2)
{
	double anglem1,anglem2;

	anglem1 = AngleO2PI(angle1+(angle2-angle1)/2.0);
	anglem2 = AngleO2PI(anglem1+PI);

	if (DistanceAngles(anglem1, angle1) < DistanceAngles(anglem2, angle1))
		return anglem1;
	else
		return anglem2;
}

inline double DistanceBtwnOrientation(const double ort1, const double ort2)
{
	double angle1[2],angle2[2],angle_min=1e10,temp;

	angle1[0] = ort1;
	angle1[1] = AngleO2PI(ort1+PI);

	angle2[0] = ort2;
	angle2[1] = AngleO2PI(ort2+PI);

	for (int i=0; i<2; i++)
		for (int j=0; j<2; j++)
		{
			temp = DistanceAngles(angle1[i],angle2[j]);
			if (temp < angle_min)
				angle_min = temp;
		}
	return angle_min;
}

inline double AverageOrientations(const double ort1, const double ort2)
{
	double angle1 = AverageAngles(ort1, ort2);
	double angle2 = AverageAngles(ort1, ort2+PI);

	if (DistanceAngles(angle1,ort1) < DistanceAngles(angle2,ort1))
		return angle1;
	else
		return angle2;
}

inline int RoundTo(const int i, const int m)
{
	return (i+2*m)%m;
}


bool GetFirstStr(int *piValue, char* lpszStr);


////////////////////////////////////////////////
////////////////////////////////////////////////
class DiscriminativeIndex
{
public:
	DiscriminativeIndex(){weight=0.0; index	= -1;};
	DiscriminativeIndex &operator=(const DiscriminativeIndex &a)
	{
		if (this != &a)
		{
			weight	= a.weight;
			index	= a.index;
		}

		return *this;
	}
	bool operator>(const DiscriminativeIndex &a){return weight>a.weight;};
	bool operator>=(const DiscriminativeIndex &a){return weight>=a.weight;};
	bool operator==(const DiscriminativeIndex &a){return weight==a.weight;};
	bool operator<(const DiscriminativeIndex &a){return weight<a.weight;};
	bool operator<=(const DiscriminativeIndex &a){return weight<=a.weight;};

public:
	double	weight;
	int		index;
};


#endif
